The present invention relates to reservoir simulation, and in particular, to methodologies for performing reservoir simulation by solving an implicit matrix equation or an implicit-IMPES matrix equation.
In an attempt to understand and predict the physical behavior of reservoirs (such as petroleum reservoirs), reservoir engineers and scientists have generated various mathematical descriptions of reservoirs and the fluids they contain. These mathematical descriptions are often expressed as coupled sets of differential equations. Since it is quite often impossible to obtain solutions of the differential equations in all but the simple cases, the differential equations are discretized in space and time, and the resulting difference equations are solved using various numerical simulation techniques. For example, the following difference equations represent the volumetric accumulation of oil and water in a particular cell (i.e. cell i) over the course of a timestep from time index n to n+1 assuming rock and fluid incompressibility in a one-dimensional reservoir:                                                                                                               (                                          λ                      o                                        )                                                        i                    +                                          1                      /                      2                                                        β                                ⁡                                  [                                                                                    (                                                  p                          o                                                )                                                                    i                        +                        1                                            α                                        -                                                                  (                                                  p                          o                                                )                                            i                      α                                                        ]                                            -                                                                    (                                          λ                      o                                        )                                                        i                    +                                          1                      /                      2                                                        β                                ⁡                                  [                                                                                    (                                                  p                          o                                                )                                            i                      α                                        -                                                                  (                                                  p                          o                                                )                                                                    i                        -                        1                                            α                                                        ]                                            +                                                (                                      q                    o                                    )                                i                β                                      =                                                            φ                  ⁢                                      xe2x80x83                                    ⁢                                      V                    i                                                                                        B                    o                                    ⁢                  Δ                  ⁢                                      xe2x80x83                                    ⁢                  t                                            ⁡                              [                                                                            (                                              S                        o                                            )                                        i                                          n                      +                      1                                                        -                                                            (                                              S                        o                                            )                                        i                    n                                                  ]                                              ,                ⁢                  
                                    (        B1        )                                                                                                      (                                      λ                    w                                    )                                                  i                  +                                      1                    /                    2                                                  β                            ⁡                              [                                                                            (                                              p                        w                                            )                                                              i                      +                      1                                        α                                    -                                                            (                                              p                        w                                            )                                        i                    α                                                  ]                                      -                                                            (                                      λ                    w                                    )                                                  i                  +                                      1                    /                    2                                                  β                            ⁡                              [                                                                            (                                              p                        w                                            )                                        i                    α                                    -                                                            (                                              p                        w                                            )                                                              i                      -                      1                                        α                                                  ]                                      +                                          (                                  q                  w                                )                            i              β                                =                                                    φ                ⁢                                  xe2x80x83                                ⁢                                  V                  i                                                                              B                  w                                ⁢                Δ                ⁢                                  xe2x80x83                                ⁢                t                                      ⁡                          [                                                                    (                                          S                      w                                        )                                    i                                      n                    +                    1                                                  -                                                      (                                          S                      w                                        )                                    i                  n                                            ]                                      ,                            (        B2        )            
where xcex94t is the timestep size;
Vi is the volume of cell i;
xcfx86 is porosity, i.e. pore volume per cell volume;
(So)i is the saturation of oil at cell i, i.e. the fraction of the pore volume occupied by oil in cell i;
(Sw)i is the saturation of water at cell i, i.e. the fraction of the pore volume occupied by water in cell i;
Bo and Bw are the formation volume factors (FVF) for oil and water respectively;
(po)ixe2x88x921, (po)i, (po)i+1 are oil pressures at cell ixe2x88x921, cell i, and cell i+1 respectively;
(pw)ixe2x88x921, (pw)i, (pw)i+1 are water pressures at cell ixe2x88x921, cell i, and cell i+1 respectively;
(qo)i is the rate of oil injection into cell i, and takes the value zero at most cells and takes a negative value at cells which interact with a depletion well;
(qw)i is the rate of water injection into cell i, and typically takes a zero value except at cells which interact with an injection or depletion well;
(x)n and (x)n+1 represent a quantity x evaluated at time indices n and n+1 respectively, where the former is known information, having been determined from previous computations, and the later is an unknown to be solved for by some computational method; and
(x)xcex1 and (x)xcex2 represent quantities which are to be evaluated at time index n or n+1 subject to user selection.
The oil transmissibility-mobility factors (xcexo)i+xc2xd and (xcexo)ixe2x88x92xc2xd are defined as                                                         (                              λ                o                            )                                      i              +                              1                /                2                                              =                                    (                              Ak                                                      x                                          i                      +                      1                                                        -                                      x                    i                                                              )                        ⁢                                          (                                  M                  o                                )                                            i                +                                  1                  /                  2                                                                    ,                            (        B3        )                                                                    (                              λ                o                            )                                      i              -                              1                /                2                                              =                                    (                              Ak                                                      x                    i                                    -                                      x                                          i                      -                      1                                                                                  )                        ⁢                                          (                                  M                  o                                )                                            i                -                                  1                  /                  2                                                                    ,                            (        B4        )            
where A is the area normal to the axis of the one-dimensional reservoir;
(Mo)i+xc2xd is the mobility of oil in transit between cell i and cell i+1;
(Mo)ixe2x88x92xc2xd is the mobility of oil in transit between cell i and cell ixe2x88x921;
xk is the position of the kth cell along the one-dimensional axis.
Similar definitions apply for the water transmissibility-mobility products (xcexw)i+xc2xd and (xcexw)ixe2x88x92xc2xd. The difference equations (B1) and (B2) above are augmented with several auxiliary relations as follows:
So+Sw=1,xe2x80x83xe2x80x83(B5)
pwxe2x88x92po=Pc(So),xe2x80x83xe2x80x83(B6)
Mo=Mo(po,So),xe2x80x83xe2x80x83(B7)
Mw=Mw(pw,Sw).xe2x80x83xe2x80x83(B8)
Relation (B5) follows from the definition of saturation. Capillary pressure Pc which is defined as the difference in pressure between water and oil is a known function of oil saturation. Oil mobility Mo is a known function of oil pressure and oil saturation. Water mobility Mw is a known function of water pressure and water saturation.
Since oil mobility Mo is a function of oil pressure po and oil saturation So, and these later variables are defined at cell centers, a question arises as to the proper means of evaluating the in-transit oil mobilities (Mo)i+xc2xd and (Mo)ixe2x88x92xc2xd. According to the midpoint weighting scheme, the in-transit oil mobility is defined to be the average of the mobilities at the two affected cells. For example,
(Mo)i+xc2xd=xc2xd(Mo)i+xc2xd(Mo)i+1,xe2x80x83xe2x80x83(B9)
where (Mo)i is evaluated using the oil saturation (So)i and oil pressure (po)i prevailing at cell i, and (Mo)i+1 is evaluated using the oil saturation (So)i+1 and oil pressure (po)i+1 prevailing at cell i+1. Alternatively, according to the upstream weighting scheme, the in transit mobility may be defined as the oil mobility at the upstream cell of the two affected cells, where the upstream cell is defined as the cell with higher pressure (since fluids flow from high pressure to low pressure). For example,                                           (                          M              o                        )                                i            +                          1              /              2                                      =                  {                                                                                                                (                                              M                        o                                            )                                        i                                    ,                                                                                                  if                    ⁢                                          xe2x80x83                                        ⁢                                                                  (                                                  p                          o                                                )                                            i                                                        ≥                                                            (                                              p                        o                                            )                                                              i                      +                      1                                                                                                                                                                                      (                                              M                        o                                            )                                                              i                      +                      1                                                        ,                                                                              otherwise                  .                                                                                        (        B10        )            
If the pressure variables and transmissibility-mobility factors in Equations (B1) and (B2) are evaluated at the new time index, i.e. xcex1=xcex2n+1, Equations (B1) and (B2) take the form                                                                                           (                                      λ                    o                                    )                                                  i                  +                                      1                    /                    2                                                                    n                  +                  1                                            ⁡                              [                                                                            (                                              p                        o                                            )                                                              i                      +                      1                                                              n                      +                      1                                                        -                                                            (                                              p                        o                                            )                                        i                                          n                      +                      1                                                                      ]                                      -                                                            (                                      λ                    o                                    )                                                  i                  -                                      1                    /                    2                                                                    n                  +                  1                                            ⁡                              [                                                                            (                                              p                        o                                            )                                        i                                          n                      +                      1                                                        -                                                            (                                              p                        o                                            )                                                              i                      -                      1                                                              n                      +                      1                                                                      ]                                      +                                          (                                  q                  o                                )                            i                              n                +                1                                              =                                                    φ                ⁢                                  xe2x80x83                                ⁢                                  V                  i                                                                              B                  o                                ⁢                Δ                ⁢                                  xe2x80x83                                ⁢                t                                      ⁡                          [                                                                    (                                          S                      o                                        )                                    i                                      n                    +                    1                                                  -                                                      (                                          S                      o                                        )                                    i                  n                                            ]                                      ,                            (        B11        )                                                                                    (                                  λ                  w                                )                                            i                +                                  1                  /                  2                                                            n                +                1                                      ⁡                          [                                                                    (                                          p                      w                                        )                                                        i                    +                    1                                                        n                    +                    1                                                  -                                                      (                                          p                      w                                        )                                    i                                      n                    +                    1                                                              ]                                -                                                    (                                  λ                  w                                )                                            i                -                                  1                  /                  2                                                            n                +                1                                      ⁡                          [                                                                    (                                          p                      w                                        )                                    i                                      n                    +                    1                                                  -                                                      (                                          p                      w                                        )                                                        i                    -                    1                                                        n                    +                    1                                                              ]                                +                                    (                              q                w                            )                        i                          n              +              1                                      =                                                            φ                ⁢                                  xe2x80x83                                ⁢                                  V                  i                                                                              B                  w                                ⁢                Δ                ⁢                                  xe2x80x83                                ⁢                t                                      ⁡                          [                                                                    (                                          S                      w                                        )                                    i                                      n                    +                    1                                                  -                                                      (                                          S                      w                                        )                                    i                  n                                            ]                                .                                    (        B12        )            
The transmissibility-mobility factors and the phase injection rates are functions of saturation and pressure, and are evaluated at the new time level n+1. Thus, Equations (B11) and (B12) are non-linear in the unknown variables
(po)ixe2x88x921n+1,(po)in+1,(po)i+1n+1,
(pw)ixe2x88x921n+1,(pw)in+1,(pw)i+1n+1,
(So)ixe2x88x921n+1 and (Sw)ixe2x88x921n+1,
(So)in+1 and (Sw)in+1,
(So)i+1n+1 and (Sw)i+1n+1.xe2x80x83xe2x80x83(B13)
Equations (B11) and (B12) may be expressed in terms of a reduced set of unknown variables using relations (B5) and (B6). For example, the variable (Sw)in+1 may be replaced by 1xe2x88x92(So)in+1. Similarly, (po)in+1 may be replaced by (po)in+1+Pc[(So)in+1]. Thus, Equations (B11) and (B12) may be expressed in terms of the following reduced set of unknown variables:
xe2x80x83(po)ixe2x88x921n+1,(po)in+1,(po)i+1n+1, (So)ixe2x88x921n+1, (So)in+1,(So)i+1n+1xe2x80x83xe2x80x83(B14)
Assuming that there are N cells in the reservoir being modeled, Equations (B11) and (B12) describe a coupled non-linear system of 2N equations (two equations per cell) with 2N unknownsxe2x80x94each cell contributes an unknown pressure (po)in+1 and an unknown saturation (So)in+1. An iterative method such as Newton""s method is generally required to solve such systems.
Let vector X be the vector of 2N unknowns for the system. Define a set of 2N functions fj, j=0, 1, 2, 3, . . . , 2Nxe2x88x921, two functions per cell, as follows. A first function f2i(X) for cell i is defined by the expression which follows from subtracting the right-hand side of Equation (B11) from the left-hand side of Equation (B11). A second function f2i+1(X) for cell i is defined by the expression which follows from subtracting the right-hand side of Equation (B12) from the left-hand side of Equation (B12). Let f: R2Nxe2x86x92R2N be the corresponding vector function whose component functions are the functions fj. The system given by Equations (B11) and (B12) may be equivalently expressed by the equation
ƒ(X)={right arrow over (0)},xe2x80x83xe2x80x83(B15)
i.e. the solution X=X* of the system given by Equations (B11) and (B12) corresponds to the zero of Equation (B15). Equation (B15) may be referred to as a fully implicit equation or a nonlinear implicit equation since none of the unknowns (B14) may be explicitly computed from known data. Thus, any method of solving equation (B15) may be referred to as a fully implicit method.
ƒ(X)={right arrow over (0)},xe2x80x83xe2x80x83(B15)
i.e. the solution X=X* of the system given by Equations (B11) and (B12) corresponds to the zero of Equation (B15). Equation (B15) may be referred to as a fully implicit equation or a nonlinear implicit equation since none of the unknowns (B14) may be explicitly computed from known data. Thus, any method of solving equation (B15) may be referred to as a fully implicit method.
Newton""s method prescribes an iterative method for obtaining the solution of Equation (B15). Given a current estimate Xk of the solution, the function ƒ is linearized in the vicinity of this current estimate:
Y=dƒ(Xk)xc2x7(Xxe2x88x92Xk)+ƒ(Xk),xe2x80x83xe2x80x83(B16)
where dƒ(Xk) represents the Jacobian matrix of the function ƒ evaluated at X=Xk, and ƒ(Xk) represents the evaluation of function ƒ at the current estimate. The next estimate Xk+1 of the solution is obtained by setting vector Y equal to zero and solving for argument X. Thus, the next estimate Xk+1 satisfies the matrix equation
dƒ(Xk)xc2x7Xk+1=dƒ(Xk)xc2x7Xkxe2x88x92ƒ(Xk).xe2x80x83xe2x80x83(B17)
By solving Equation (B17) for successively increasing values of the index k, a sequence of estimates Xo, X1, X2, . . . , Xk, . . . is obtained which converge to the solution of the nonlinear system (B15).
Equation (B17) is referred to herein as an implicit matrix equation. A linear equation solver is used to solve the implicit matrix equation (B17). The right-hand side vector dƒ(Xk)xc2x7Xkxe2x88x92ƒ(Xk) and the Jacobian matrix dƒ(Xk) are supplied as input data to the linear solver. The linear solver returns the solution vector Xk+1 of the implicit matrix equation (B17). The computational effort of a Newton""s method solution of the nonlinear implicit equation (B15) depends on (a) the number of Newton iterations to achieve convergence of the sequence Xk, (b) the average computational effort expended by the linear solver to solve the implicit matrix equation (B17), and (c) the computational effort required to update the matrix equation as improved solutions are obtained. While most of the computational effort per Newton iteration is associated with solving the matrix equation, the effort required to update the matrix equation is also significant. Thus, any improvement in the computational efficiency of the linear equation solver will have a corresponding effect on the efficiency of the Newton""s method solution.
As described above, the nonlinear implicit equation (B15) arises from the choices xcex1=xcex2n+1 in Equations (B1) and (B2) above. Another plausible set of choices is given by xcex1=n+1 and xcex2=n , whereupon Equations (B1) and (B2) take the form                                                                                           (                                      λ                    o                                    )                                                  i                  +                                      1                    /                    2                                                  n                            ⁡                              [                                                                            (                                              p                        o                                            )                                                              i                      +                      1                                                              n                      +                      1                                                        -                                                            (                                              p                        o                                            )                                        i                                          n                      +                      1                                                                      ]                                      -                                                            (                                      λ                    o                                    )                                                  i                  -                                      1                    /                    2                                                  n                            ⁡                              [                                                                            (                                              p                        o                                            )                                        i                                          n                      +                      1                                                        -                                                            (                                              p                        o                                            )                                                              i                      -                      1                                                              n                      +                      1                                                                      ]                                      +                                          (                                  q                  o                                )                            i              n                                =                                                    φ                ⁢                                  xe2x80x83                                ⁢                                  V                  i                                                                              B                  o                                ⁢                Δ                ⁢                                  xe2x80x83                                ⁢                t                                      ⁡                          [                                                                    (                                          S                      o                                        )                                    i                                      n                    +                    1                                                  -                                                      (                                          S                      o                                        )                                    i                  n                                            ]                                      ,                            (        B18        )                                                                                    (                                  λ                  w                                )                                            i                +                                  1                  /                  2                                            n                        ⁡                          [                                                                    (                                          p                      w                                        )                                                        i                    +                    1                                                        n                    +                    1                                                  -                                                      (                                          p                      w                                        )                                    i                                      n                    +                    1                                                              ]                                -                                                    (                                  λ                  w                                )                                            i                -                                  1                  /                  2                                            n                        ⁡                          [                                                                    (                                          p                      w                                        )                                    i                                      n                    +                    1                                                  -                                                      (                                          p                      w                                        )                                                        i                    -                    1                                                        n                    +                    1                                                              ]                                +                                    (                              q                w                            )                        i            n                          =                                                            φ                ⁢                                  xe2x80x83                                ⁢                                  V                  i                                                                              B                  w                                ⁢                Δ                ⁢                                  xe2x80x83                                ⁢                t                                      ⁡                          [                                                                    (                                          S                      w                                        )                                    i                                      n                    +                    1                                                  -                                                      (                                          S                      w                                        )                                    i                  n                                            ]                                .                                    (        B19        )            
The saturations and pressures at time-index n comprise known data (having been determined from previous computations). Thus, the transmissibility-mobility functions evaluated at time-index n comprise known constants. Equations (B18) and (B19) are therefore linear in the unknown variables
(po)ixe2x88x921n+1, (po)in+1, (po)i+1n+1,
(pw)ixe2x88x921n+1, (pw)in+1, (pw)i+1n+1,
(So)in+1 and (Sw)in+1.xe2x80x83xe2x80x83(B20)
One method for solving the linear system of Equations (B18) and (B19), i.e. the so called Implicit-Pressure Explicit-Saturation (IMPES) method, is motivated by the following reduction of Equations (B18) and (B19). Since the saturation variables obey relation (B5), the Equations (B18) and (B19) may be combined so as to eliminate the unknown saturation variables. In particular, Equation (B18) may be multiplied by the oil formation volume factor Bo, and Equation (B19) may be multiplied by the water formation volume factor Bw. The resulting equations may be added together to generate the following linear equation involving only the pressure unknowns:
Bo(xcexo)i+xc2xdn[(po)i+1n+1xe2x88x92(po)in+1]xe2x88x92Bo(xcexo)ixe2x88x92xc2xdn[(po)in+1xe2x88x92(po)ixe2x88x921n+1]
+Bw(xcexw)i+xc2xdn[(pw)i+1n+1xe2x88x92(pw)in+1]xe2x88x92Bw(xcexw)ixe2x88x92xc2xdn[(pw)in+1xe2x88x92(pw)ixe2x88x921n+1]
+Bo(qo)in+Bw(qw)in=0.xe2x80x83xe2x80x83(B21)
Equation (B21) is referred to herein as an IMPES pressure equation. The capillary pressure relation (B6) may be used to eliminate the water pressure unknowns under the assumption that capillary pressure does not change during the timestep:
(pw)jn+1=(po)jn+1+Pc[(So)jn],xe2x80x83xe2x80x83(B22)
where j represents an arbitrary cell index. When Equation (B21) is written for all N cells in the reservoir, the ensuing system, herein referred to as the IMPES pressure system, has N equations and N unknownsxe2x80x94one unknown pressure (po)jn+1 per cell. Because the IMPES pressure system is linear and has fewer equations and unknowns it may be solved much faster than the fully implicit system (B15).
Again a linear equation solver may be invoked to solve the IMPES pressure system. The solution vector pn+1 of the IMPES pressure system specifies the pressure (po)in+1 at every cell in the reservoir. The unknown saturations (So)in+1 and (Sw)in+1 in Equations (B18) and (B19) may be determined by substituting the pressure solution values (po)jn+1 into the left-hand sides of Equations (B18) and (B19). Since the saturations (So)in and (Sw)in are known from previous computations, the unknown saturations (So)in+1 and (Sw)in+1 may be computed explicitly. Thus, the IMPES procedure involves two steps: a first step in which pressures are computed implicitly as the solution of a linear system; and a second step in which saturations are computed explicitly based on the pressure solution.
The example of a one-dimensional model discussed above represents a greatly simplified description of a complicated physical situation. More realistic models involve (a) a two-dimensional or three-dimensional array of cells, (b) more than two conserved species, (c) more than two phases, (d) compressible fluids and/or rock substrate, (e) non-uniform cell geometry and spacing, etc. In addition, the difference equations of the reservoir model may not necessarily arise from a fluid volume balance. In other approaches, difference equations may be obtained by performing, e.g., mass or energy balances. While pressure is quite often one of the variables being solved for at each cell, the remaining variables need not necessarily be saturations. For example, in other formulations, the remaining variables may be mole fractions, masses, or other quantities.
Given a reservoir with M conserved species, a conservation law may be invoked to write a set of M difference equations describing the physical behavior of each of the conserved species at a generic cell i. (The use of a single index i to denote a generic cell does not necessarily imply that the reservoir model is one-dimensional.) The set of equations may generally be expressed in terms of the pressure Pb of some base species (often oil), and (Mxe2x88x921) complementary variables such as saturations, mole fractions, masses, etc. These complementary variables will be referred to herein as generalized saturations.
The discussion of the fully implicit method and the IMPES method presented above generalizes to more realistic models. The M difference equations for the generic cell i generally include functions such as mobility, formation volume factor, pore volume, injection rate etc., which depend on pressure and/or the generalized saturations (i.e. complementary variables). The fully implicit equations result from evaluating such functions at the new time index n+1. The fully implicit equations are generally non-linear, and thus, require an iterative method such as Newton""s method for their solution.
The IMPES formulation starts from evaluating functions of pressure and/or the complementary variables at the old time index n. Thus, the M difference equations particularize to a set of linear equations in the unknown pressures and unknown generalized saturations. An auxiliary relation analogous to relation (B5) may be used to combine the set of linear equations into a single equation which involves only the pressure unknowns. This single equation is commonly referred to as the IMPES pressure equation. The IMPES pressure equation may be solved by calling a linear equation solver. The pressure solution is then substituted into the original set of linear equations, and the generalized saturations are computed explicitly.
Both the fully implicit method and the IMPES method aim at generating values for the base pressure and the generalized saturations at the new time index n+1 for each cell in the reservoir. However, because the IMPES method is less stable than the fully implicit method (FIM), the timestep xcex94tIMPES used in the IMPES method is generally significantly smaller than the timestep xcex94tFIM used in the fully implicit method. While the single-timestep computational effort CEIMPES of IMPES is much smaller than the single-timestep computational effort CEFIM of the fully implicit method, it is quite often the case that the ratio       Δ    ⁢          xe2x80x83        ⁢          t      FIM            Δ    ⁢          xe2x80x83        ⁢          t      IMPES      
of timestep sizes is larger than the ratio       CE    FIM        CE    IMPES  
of computational efforts. Thus, any advantage gained by the single-timestep efficiency of the IMPES method is counteracted by the necessity of performing a large number of IMPES timesteps to cover a timestep of the fully implicit method.
The IMPES method is one method in a general class of methods commonly referred to as sequential methods. A sequential method involves a two-step procedure: a first step in which unknown pressures are determined, and a second step in which comlementary unknowns (i.e. unknowns other than pressure) are determined using the pressure solution obtained in the first step.
Another sequential method, commonly referred to as the total velocity sequential semi-implicit (TVSSI) method has received significant use since it was originally developed by Spillette et al. circa 1970. The TVSSI method is described in the following paper by Spillette, A. G., Hillestad, J. G., and Stone, H. L.: xe2x80x9cA High-Stability Sequential Solution Approach to Reservoir Simulation,xe2x80x9d SPE 4542 presented at the 1973 SPE Annual Meeting, Las Vegas, September 30-October 3. This paper is hereby incorporated by reference.
Similar to the IMPES method, the TVSSI method has the advantage of reduced computational effort per timestep as compared to the fully implicit method. However, the TVSSI method is far more stable than the IMPES method. The increased stability implies that the timestep xcex94tTVSSI of the TVSSI method may be significantly larger than the IMPES timestep xcex94tIMPES. The TVSSI method comprises two major steps: (i) solving the IMPES pressure system; and (ii) solving a set of coupled saturation equations for the generalized saturations. Since the IMPES pressure equation involves a single unknown (i.e. pressure) at each cell, step (i) requires significantly less work than solving the set of fully implicit equations. In addition, since the set of coupled saturation equations does not have the elliptic nature of the IMPES pressure equation or the set of fully implicit equations, the saturation solution converges rapidly. Overall, the single-timestep computational effort CETVSSI for the TVSSI method is typically a half to a fifth that of the fully implicit computations.
The TVSSI method is not as stable as the fully implicit method. In some problems, the ratio       Δ    ⁢          xe2x80x83        ⁢          t      FIM            Δ    ⁢          xe2x80x83        ⁢          t      TVSSI      
of timestep sizes is larger than the ratio       CE    FIM        CE    TVSSI  
of computational efforts. In other words, the single timestep computational efficiency of the TVSSI method relative to the fully implicit method is more than offset by the necessity of performing multiple timesteps of the TVSSI method to cover a timestep of the fully implicit method.
Overall, the fully implicit method seems to be more desirable than the TVSSI method, in part because it is more trouble-free. However, the total velocity equations contain a certain power that enables the success, albeit not universal, of the TVSSI method. This power has yet to be fully appreciated and harnessed. Thus, there exists a need for a reservoir simulation method which may more effectively capture this power inherent in the total velocity equations.
One prior-art method used to lower the cost of reservoir simulations is the so called adaptive implicit method (AIM). The adaptive implicit method is based on the recognition that the implicit formulation is required at only a fraction of the cells in the reservoir model. If the implicit formulation can be applied only where it is needed, with the IMPES formulation being used at the remaining cells, significant reductions in computational effort may be obtained. The adaptive implicit method determines dynamically which cells require implicit formulation. As the simulation progresses in time, a particular cell may switch back and forth between IMPES formulation and implicit formulation.
In a related prior-art method, referred to as static variable implicitness, the assignment of IMPES or implicit formulation to each cell in the reservoir remains fixed through the simulation.
Although the adaptive implicit method and variable implicit method are computationally more efficient than the fully implicit method, they are still significantly time consuming. Thus, there exists a need for improved methods for performing adaptive and variable implicit reservoir simulations.
The present invention comprises a method for performing reservoir simulation by solving a mixed implicit-IMPES matrix (MIIM) equation. The MIIM equation arises from a Newton iteration of a variable implicit reservoir model. The variable implicit reservoir model comprises a plurality of cells including both implicit cells and IMPES cells. The MIIM equation includes a scalar IMPES equation for each of the IMPES cells and a set of implicit equations for each of the implicit cells.
In one embodiment, the method for performing reservoir simulation comprises: (a) constructing a global IMPES pressure matrix equation from the MIIM equation; (b) determining coefficients for a set of saturation equations at the implicit cells by using a total velocity constraint at the implicit cells; (c) solving the global IMPES pressure matrix equation for pressure changes; (d) computing first residuals at the implicit cells in response to the pressure changes; (e) solving the set of saturation equations (formed from the coefficients and first residuals) for saturation changes at the implicit cells; (f) computing second residuals at the implicit cells and at a subset of the IMPES cells that are in flow communication with any of the implicit cells in response to the saturation changes. Steps (b) through (f) may be repeated until the second residuals satisfy a convergence condition. A final solution estimate may be computed for the MIIM equation from the pressures changes and the saturation changes after the convergence condition is satisfied. The final solution estimate may be used by a reservoir simulator to determine behavior of the reservoir model at a future discrete time value.
The global IMPES pressure matrix equation may be constructed from the MIIM equation by (i) manipulating the set of implicit equations at each implicit cell to generate a corresponding IMPES pressure equation, and (ii) concatenating the IMPES pressure equations for the IMPES cells and the IMPES pressure equations for the implicit cells. Note the IMPES pressure equations for the IMPES cells are provided by the MIIM equations.
In a second embodiment, the method for performing reservoir simulation comprises: (a) constructing a global IMPES pressure equation from the MIIM equation; (b) solving the global IMPES pressure equation for pressure changes; (c) computing first residuals at the implicit cells in response to the pressure changes; (d) determining improved saturations and improved pressures by performing one or more iterations with a selected preconditioner at the implicit cells; and (e) computing second residuals at the implicit cells and at a subset of the IMPES cells that are in flow communication with any of the implicit cells in response to the improved saturations and improved pressures. Steps (b) through (e) may be repeated until a convergence condition based on the second residuals is satisfied. A final solution estimate for the MIIM equation may be computed from the pressure changes, improved saturations and improved pressures after the convergence condition is satisfied. The final solution estimate may be used to determine behavior of the reservoir model at a future discrete time value.
In a third embodiment, the method for performing reservoir simulation comprises: (a) constructing a global IMPES pressure equation from the MIIM equation; (b) solving the global IMPES pressure equation for pressure changes; (c) computing first residuals at the implicit cells in response to the pressure changes; (d) solving an implicit system comprising the set of implicit equations associated with each of the implicit cells for improved saturations and improved pressures at the implicit cells using the first residuals at the implicit cells; and (e) computing second residuals for a subset of the IMPES cells which are in flow communication with any of the implicit cells. Steps (b) through (e) may be iterated until a convergence condition is satisfied based on the second residuals. The final solution estimate for the MIIM equation may be computed based on the improved saturations and improved pressures after the convergence condition is satisfied. In solving the implicit system, cell pressures for fringe IMPES cells (i.e. the IMPES cells which are in flow communication with any implicit cell) are held fixed at those values determined in the pressure solution of step (b).